package variance;

import java.lang.Math;
import yieldCurve.*;

public class BlackScholes {
	// S, expiry, K, P are read from benchmark European options
	YieldCurve yieldCurve;
	double S; // Stock price now
	Date expiry; // Expiry date
	Date today; // today
	double T; // delta T from today to expiry date
	double K; // strike price
	double P; // premium
	double r; // continuous compound rate at expiry date
	
	// Private variables used for Black-Scholes formula
	static final double beta = 0.2316419;
	static final double a1 = 0.319381530;
	static final double a2 = -0.356563782;
	static final double a3 = 1.781477937;
	static final double a4 = -1.821255978;
	static final double a5 = 1.330274429;
	
	// initial guess of sigma x0;
	static double initialGuess = 0.1;
	
	/** Class constructor
	 * 
	 * @param curve A populated YieldCurve
	 * @param S current underlying price
	 * @param today Start date
	 * @param expiry Expiry date
	 * @param K strike price
	 * @param P benchmark option price
	 */
	public BlackScholes(YieldCurve curve, double S, Date today, Date expiry, double K, double P) {
		yieldCurve = curve;
		this.S = S;
		this.expiry = expiry;
		this.today = Date.businessDay(today);
		this.K = K;
		this.P = P;
		r = yieldCurve.getContinuousCompound(expiry);
		T = Date.getDeltaT(today, expiry);
		
	}
	
	/**
	 * The Approximate function of the standard normal distribution N(x)
	 * @param x
	 * @return probability of normal distribution given x
	 */
	private double N(double x) {
		double n_dash;
		double k = 1 / (1 + beta * x);
		
		n_dash = 1 / Math.sqrt(2 * Math.PI) * Math.exp(- x * x /2);
		if (x >= 0) {
			return 1 - n_dash * (a1*k + a2*Math.pow(k, 2) + a3*Math.pow(k, 3) + a4*Math.pow(k, 4) + a5*Math.pow(k, 5));
		} else {
			return 1 - N(-x);
		}
	}
	/**
	 * The object equation for a European call
	 * f = C - premium
	 * @param sigma
	 * @return f for Newton's method
	 */
	public double f(double sigma) {
		double d1 = (Math.log(S/K) + (r + sigma*sigma/2) * T) / (sigma * Math.sqrt(T));
		double d2 = d1 - sigma*Math.sqrt(T);
		double C = S * N(d1) - K * Math.exp(-r*T) * N(d2);
		return C - P;
	}
	
	/**
	 * Get volatility/sigma using Newton's method
	 * @return volatility
	 */
	public double getSigma() {
		// x0 is initial guess of sigma
		// To make d1 and d2 reusable in f and f', write them here. x0, x1 are just
		// the same as sigma in Black-Schole formula
		double x0 = initialGuess;
		double x = x0;
		double f;
		double f_dash;
		double d1;
		double d2;
		double C;
		
		double epsilon = 0.0001;
		
		do {
			x0 = x;
			// get f
			d1 = (Math.log(S/K) + (r + x0*x0/2) * T) / (x0 * Math.sqrt(T));
			d2 = d1 - x0*Math.sqrt(T);
			C = S * N(d1) - K * Math.exp(-r*T) * N(d2);
			f = C - P;
			
			//get f's first order derivative f_dash
			// Approximate method given in the lecture notes
			f_dash = S * Math.sqrt(T) * (1/Math.sqrt(2*Math.PI)) * Math.exp(-d1*d1/2);

/*			My precise method...
			double dNdd1;
			double dd1dsigma;
			double dNdd2;
			double dd2dsigma;
			dNdd1 = 1 / Math.sqrt(2*Math.PI) * Math.exp(-d1*d1/2);
			dd1dsigma = (x0*x0*T/2 - Math.log(S/K) - r*T)/ (x0*x0*Math.sqrt(T));
			dNdd2 = 1 / Math.sqrt(2*Math.PI) * Math.exp(-d2*d2/2);
			dd2dsigma = dd1dsigma - Math.sqrt(T);
			f_dash = S * dNdd1 * dd1dsigma - K*Math.exp(-r*T)*dNdd2*dd2dsigma;
*/			

			// Newton's method
			x = x0 - f / f_dash;
		} while (Math.abs((x - x0)/ x) > epsilon);
		return x;
	}

	/** Test method
	 * 
	 * @param args Parameters to executable (not required or supported)
	 */
	public static void main(String[] args) {
		double sigma;
		Date today = new Date(2009,5,1);
		Date expiry = new Date(2011, 1, 22);
		YieldCurve yieldCurve = new YieldCurve("curveDataInput.csv", "curveSpec.csv", 2, today);
		BlackScholes bc = new BlackScholes(yieldCurve, 89.31, today, expiry, 95, 25);
		sigma = bc.getSigma();
		System.out.print("sigma=" + sigma + "\n");
		System.out.print("N(sigma)=" + bc.N(sigma) + "\n");
		System.out.print("f(sigma)=" + bc.f(sigma)+ "\n");
	}
}
